#include "Plotter.h"
#include "TMath.h"
#include "TFile.h"
#include "TAxis.h"

TCanvas* plot_mtop(TString hist_fpath, const int ijet, const int ierr) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "mtop_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    TFile *f = new TFile(hist_fpath);

    TH1D * h_data    = (TH1D*) f->Get("data10_mtop"     + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mtop"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mtop"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mtop"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mtop"   + ht_ext); h_zjets->SetDirectory(0);

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mtop_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mtop_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */
    TString prefix_hqcd_name = "data10_mtop_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("MeV");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("MeV");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    // adding up histogram for making stack plot
    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    h_ttbar->Add(h_qcd);

    // deciding the first drawing histogram (data or MC?)
    Double_t max_data = 0;
    Double_t max_mc   = 0;

    for ( int i=1; i<=h_data->GetNbinsX(); i++) {
        Double_t cbin = h_data->GetBinContent(i);
        if ( cbin > max_data ) max_data = cbin; 
    }

    for ( int i=1; i<=h_ttbar->GetNbinsX(); i++) {
        Double_t cbin = h_ttbar->GetBinContent(i);
        if ( cbin > max_mc ) max_mc = cbin;
    }

    if ( max_data > max_mc ) {
        h_data->Draw("E1");
        h_ttbar->Draw("HIST SAME");
    } else {
        h_ttbar->Draw("HIST");
    }

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    h_qcd->Draw("HIST SAME");

    h_qcd_err->SetMarkerSize(0);
    h_qcd_err->SetFillStyle(3001);
    h_qcd_err->SetFillColor(kBlack);
    h_qcd_err->Draw("E2 SAME");

    h_data->Draw("E1 SAME");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    legend->AddEntry(h_qcd,"QCD","f");
    legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    cout << "KS test: " << KSTest << endl;

    TString ks_str = "KS: ";
    ks_str += KSTest;

    myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    //f->Close();

    return c;
}

TCanvas* plot_eff_l2t_data(TString hist_fpath, const int ijet) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "eff_l2t_data_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    TFile *f = new TFile(hist_fpath);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    TH1D * h_tight = (TH1D*) f->Get("data10_d0sig_tight" + ht_ext); h_tight->SetDirectory(0);
    TH1D * h_loose = (TH1D*) f->Get("data10_d0sig_loose" + ht_ext); h_loose->SetDirectory(0);
 
    TGraphAsymmErrors *g = EfficiencyPlot(h_loose, h_tight);
    g->SetMarkerColor(kBlack);
    g->SetMarkerStyle(20);
    g->Draw("P");

    g->GetXaxis()->SetTitle("d0 significance");
    g->GetYaxis()->SetTitle("efficiency");
    g->GetXaxis()->SetLabelSize(.03);
    g->GetYaxis()->SetLabelSize(.03);
    g->GetXaxis()->SetTitleSize(.03);
    g->GetYaxis()->SetTitleSize(.03);
    g->GetXaxis()->SetTitleOffset(1.5);
    g->GetYaxis()->SetTitleOffset(1.5);

    g->SetLineColor(kBlue);
    g->SetMarkerColor(kBlue);
    g->SetMarkerStyle(22);
    g->Draw("ALP");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(g,"data (36.6 pb^{-1})","p");

    legend->Draw();

    //myText2(0.6, 0.8, 1, .03, plot_title);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}

TCanvas* plot_eff_l2t_mc(TString hist_fpath, const int ijet, const bool inc_mcqcd, const bool inc_data) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "eff_l2t_mc_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    TFile *f = new TFile(hist_fpath);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    TH1D *h_tight_mix = (TH1D*) f->Get("mc_mix_d0sig_tight" + ht_ext); h_tight_mix->SetDirectory(0);
    TH1D *h_loose_mix = (TH1D*) f->Get("mc_mix_d0sig_loose" + ht_ext); h_loose_mix->SetDirectory(0);
 
    TGraphAsymmErrors *g_mix = EfficiencyPlot(h_loose_mix, h_tight_mix);

    g_mix->GetXaxis()->SetTitle("d0 significance");
    g_mix->GetYaxis()->SetTitle("efficiency");
    g_mix->GetXaxis()->SetLabelSize(.03);
    g_mix->GetYaxis()->SetLabelSize(.03);
    g_mix->GetXaxis()->SetTitleSize(.03);
    g_mix->GetYaxis()->SetTitleSize(.03);
    g_mix->GetXaxis()->SetTitleOffset(1.5);
    g_mix->GetYaxis()->SetTitleOffset(1.5);

    g_mix->SetLineColor(kBlue);
    g_mix->SetMarkerColor(kBlue);
    g_mix->SetMarkerStyle(22);
    g_mix->Draw("ALP");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);

    legend->AddEntry(g_mix,"MC mix","lp");

    if ( inc_mcqcd ) {

        TH1D *h_tight_qcd = (TH1D*) f->Get("mc_qcd_d0sig_tight" + ht_ext); h_tight_qcd->SetDirectory(0);
        TH1D *h_loose_qcd = (TH1D*) f->Get("mc_qcd_d0sig_loose" + ht_ext); h_loose_qcd->SetDirectory(0);
 
        TH1D *h_tight_qcdbb = (TH1D*) f->Get("mc_qcdbb_d0sig_tight" + ht_ext); h_tight_qcdbb->SetDirectory(0);
        TH1D *h_loose_qcdbb = (TH1D*) f->Get("mc_qcdbb_d0sig_loose" + ht_ext); h_loose_qcdbb->SetDirectory(0);

        h_loose_qcd->Add(h_loose_qcdbb);
        h_tight_qcd->Add(h_tight_qcdbb);

        TGraphAsymmErrors *g_qcd = EfficiencyPlot(h_loose_qcd, h_tight_qcd);

        g_qcd->SetLineColor(kRed);
        g_qcd->SetMarkerColor(kRed);
        g_qcd->SetMarkerStyle(21);
        g_qcd->Draw("LP");

        legend->AddEntry(g_qcd,"MC QCD + QCDbb (filtmu)","lp");
    }

    if ( inc_data ) {
        TH1D * h_tight = (TH1D*) f->Get("data10_d0sig_tight" + ht_ext); h_tight->SetDirectory(0);
        TH1D * h_loose = (TH1D*) f->Get("data10_d0sig_loose" + ht_ext); h_loose->SetDirectory(0);
 
        TGraphAsymmErrors *g = EfficiencyPlot(h_loose, h_tight);
        g->SetMarkerColor(kBlack);
        g->SetMarkerStyle(20);
        g->Draw("P");

        legend->AddEntry(g,"data (36.6 pb^{-1})","p");
    }
    legend->Draw();

    //myText2(0.6, 0.8, 1, .03, plot_title);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}


TCanvas* plot_nbjets(TString hist_fpath, const int ijet, const bool logy, const bool qcd, const int ierr) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "nbjets25_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    TFile *f = new TFile(hist_fpath);

    TH1D * h_data    = (TH1D*) f->Get("data10_nbjets25"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_nbjets25"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_nbjets25"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_nbjets25"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_nbjets25"   + ht_ext); h_zjets->SetDirectory(0);

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_nbjets25_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_nbjets25_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */
    TString prefix_hqcd_name = "data10_nbjets25_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }


    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("NBjets");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("NBjets");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();
    if (logy) c->SetLogy();

    // adding up histogram for making stack plot
    h_data->Draw("E1");

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd) h_ttbar->Add(h_qcd);
    h_ttbar->Draw("HIST SAME");

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    //Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    //cout << "KS test: " << KSTest << endl;

    //TString ks_str = "KS: ";
    //ks_str += KSTest;

    //myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    // showing each bin value data:qcd
    for (Int_t i=0; i < h_data->GetNbinsX() + 1; i++) {
        Double_t N_sel     = h_data->GetBinContent(i);
        Double_t N_qcd     = h_qcd->GetBinContent(i);
        Double_t N_qcd_err = h_qcd->GetBinError(i);
        cout << i << ":" << N_sel << " " << N_qcd << " +- " << N_qcd_err << endl;
    }

    delete f;

    return c;
}

TCanvas* plot_njets(TString hist_fpath, const int ijet, const bool logy, const bool qcd, const int ierr) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "njets25_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    TFile *f = new TFile(hist_fpath);

    TH1D * h_data    = (TH1D*) f->Get("data10_njets"       + ht_ext); h_data->SetDirectory(0); set_jet_binlabel(h_data);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_njets"   + ht_ext); h_ttbar->SetDirectory(0);set_jet_binlabel(h_ttbar);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_njets"    + ht_ext); h_stop->SetDirectory(0); set_jet_binlabel(h_stop);
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_njets"   + ht_ext); h_wjets->SetDirectory(0);set_jet_binlabel(h_wjets);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_njets"   + ht_ext); h_zjets->SetDirectory(0);set_jet_binlabel(h_zjets);

    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();set_jet_binlabel(h_qcd);

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_njets_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_njets_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */

    TString prefix_hqcd_name = "data10_njets_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0); set_jet_binlabel(h_qcd_err);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("Njets");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("Njets");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();
    if (logy) c->SetLogy();

    // adding up histogram for making stack plot
    h_data->Draw("E1");

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd) h_ttbar->Add(h_qcd);
    h_ttbar->Draw("HIST SAME");

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    //Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    //cout << "KS test: " << KSTest << endl;

    //TString ks_str = "KS: ";
    //ks_str += KSTest;

    //myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    // showing each bin value data:qcd
    for (Int_t i=0; i < h_data->GetNbinsX() + 1; i++) {
        Double_t N_sel     = h_data->GetBinContent(i);
        Double_t N_qcd     = h_qcd->GetBinContent(i);
        Double_t N_qcd_err = h_qcd->GetBinError(i);
        cout << i << ":" << N_sel << " " << N_qcd << " +- " << N_qcd_err << endl;
    }

    delete f;

    return c;
}

TCanvas* plot_mwt(TString hist_fpath, const int ijet, const bool qcd, const bool mc_qcd, const double mc_qcd_scale, const int ierr) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "mwt_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    ostringstream st_scale_mc2d;
    st_scale_mc2d << mc_qcd_scale;

    TFile *f = new TFile(hist_fpath);

    TH1D * h_data    = (TH1D*) f->Get("data10_mwt"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mwt"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mwt"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mwt"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mwt"   + ht_ext); h_zjets->SetDirectory(0);

    TH1D * h_mcqcd   = h_data->Clone(); h_mcqcd->Reset();
    TH1D * h_mcqcdbb = h_data->Clone(); h_mcqcdbb->Reset();

    if ( mc_qcd ) {
        h_mcqcd   = (TH1D*) f->Get("mc_qcd_mwt"     + ht_ext); h_mcqcd->SetDirectory(0);
        h_mcqcdbb = (TH1D*) f->Get("mc_qcdbb_mwt"   + ht_ext); h_mcqcdbb->SetDirectory(0);

        h_mcqcd->Scale(mc_qcd_scale);
        h_mcqcdbb->Scale(mc_qcd_scale);
    }

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mwt_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mwt_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */

    TString prefix_hqcd_name = "data10_mwt_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("MeV");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("MeV");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);
    h_qcd_err->SetFillColor(kOrange+4);

    h_mcqcd->SetFillColor(kOrange+2);
    h_mcqcdbb->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    // adding up histogram for making stack plot

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd)    h_ttbar->Add(h_qcd);
    if (mc_qcd) {
        h_ttbar->Add(h_mcqcd);
        h_ttbar->Add(h_mcqcdbb);
    }
    // deciding the first drawing histogram (data or MC?)
    Double_t max_data = 0;
    Double_t max_mc   = 0;

    for ( int i=1; i<=h_data->GetNbinsX(); i++) {
        Double_t cbin = h_data->GetBinContent(i);
        if ( cbin > max_data ) max_data = cbin; 
    }

    for ( int i=1; i<=h_ttbar->GetNbinsX(); i++) {
        Double_t cbin = h_ttbar->GetBinContent(i);
        if ( cbin > max_mc ) max_mc = cbin;
    }

    if ( max_data > max_mc ) {
        h_data->Draw("E1");
        h_ttbar->Draw("HIST SAME");
    } else {
        h_ttbar->Draw("HIST");
    }

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    if (mc_qcd) {
        h_stop->Add(h_mcqcd);
        h_stop->Add(h_mcqcdbb);
    }
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    if (mc_qcd) {
        h_wjets->Add(h_mcqcd);
        h_wjets->Add(h_mcqcdbb);
    }
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    if (mc_qcd) {
        h_zjets->Add(h_mcqcd);
        h_zjets->Add(h_mcqcdbb);
    }
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    if (mc_qcd) {
        h_mcqcd->Draw("HIST SAME");
        h_mcqcdbb->Draw("HIST SAME");
    }

    h_data->Draw("E1 SAME");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    if (mc_qcd) {
        legend->AddEntry(h_mcqcd,"QCD*" + TString(st_scale_mc2d.str()),"f");
        legend->AddEntry(h_mcqcdbb,"QCDbb*" + TString(st_scale_mc2d.str()),"f");
    }
    legend->Draw();

    // KS test
    Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    cout << "KS test: " << KSTest << endl;

    TString ks_str = "KS: ";
    ks_str += KSTest;

    myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}

TCanvas * cal_perf_syserr(TString hist_fpath, const int ijet, const int ierr, const bool logy) {

    TString ht_ext = "_j";
    ht_ext += ijet;

    ostringstream stm_title;
    stm_title << "sys_err_perf" << ht_ext;
    if ( hist_fpath.Contains("_btag") ) stm_title << "_btag";
    if ( hist_fpath.Contains("_tria") ) stm_title << "_tria";
    TString plot_title = TString(stm_title.str());
    
    TFile * f = new TFile( hist_fpath );

    TString ht_mid = "";
    switch(ierr) {
        case 1:
            ht_mid += "_estat";
            break;
      
        case 2:
            ht_mid += "_esyst";
            break;

        default:
            break;
    }

    ht_mid += ht_ext;

    TH1D  * h_mc_qcd     = (TH1D*) f->Get("mc_qcd_njets"     + ht_mid); h_mc_qcd->SetDirectory(0);
    TH1D  * h_mc_qcdbb   = (TH1D*) f->Get("mc_qcdbb_njets"   + ht_mid); h_mc_qcdbb->SetDirectory(0);
    TH1D  * h_mc_qcd_est = (TH1D*) f->Get("mc_mix_njets_qcd" + ht_mid); h_mc_qcd_est->SetDirectory(0);

    h_mc_qcd->Add(h_mc_qcdbb);
   
    vector<Double_t>  qcd_mc;          qcd_mc.resize(10);
    vector<Double_t>  qcd_mc_err;      qcd_mc_err.resize(10);
    vector<Double_t>  qcd_mc_est;      qcd_mc_est.resize(10);
    vector<Double_t>  qcd_mc_est_err;  qcd_mc_est_err.resize(10);

    // getting full spectrum (0-9 jets) of the jet multiplicity
    int s = (int) qcd_mc.size();

    qcd_est_by_jets(h_mc_qcd,     &qcd_mc,     &qcd_mc_err);
    qcd_est_by_jets(h_mc_qcd_est, &qcd_mc_est, &qcd_mc_est_err);

    for (int i=0; i < s; i++) {
        if ( qcd_mc.at(i) > 0 ) {
            Double_t r = qcd_mc.at(i) / qcd_mc_est.at(i);
            cout << "ijet=" << i << " " << qcd_mc.at(i) << " " << qcd_mc_est.at(i) << " " << r << endl;
        } else {
            cout << "ijet=" << i << " " << qcd_mc.at(i) << " " << qcd_mc_est.at(i) << endl;
        }
    }

    // making histogram up to 4 jet inclusive bin (> 4 jets are summed up into 4 jet inclusive bin)
    qcd_mc.resize(5);
    qcd_mc_err.resize(5);
    qcd_mc_est.resize(5);
    qcd_mc_est_err.resize(5);

    qcd_est_by_jets(h_mc_qcd,     &qcd_mc,     &qcd_mc_err);
    qcd_est_by_jets(h_mc_qcd_est, &qcd_mc_est, &qcd_mc_est_err);

    s = (int) qcd_mc.size();

    TH1D * hp_mc_qcd = new TH1D("hp_mc_qcd", "hp_mc_qcd", 5, 0, 5); hp_mc_qcd->SetDirectory(0);
    TH1D * hp_mm_qcd = new TH1D("hp_mm_qcd", "hp_mc_qcd", 5, 0, 5); hp_mc_qcd->SetDirectory(0);

    for (int i=0; i < s; i++) {
        if ( qcd_mc.at(i) > 0 ) {
            Double_t r = qcd_mc.at(i) / qcd_mc_est.at(i);
            cout << "ijet=" << i << " " << qcd_mc.at(i) << " " << qcd_mc_est.at(i) << " " << r << endl;
        } else {
            cout << "ijet=" << i << " " << qcd_mc.at(i) << " " << qcd_mc_est.at(i) << endl;
        }

        if ( i > 0 ) {
            hp_mc_qcd->SetBinContent(i+1, qcd_mc.at(i));
            hp_mc_qcd->SetBinError(i+1, qcd_mc_err.at(i));
            hp_mm_qcd->SetBinContent(i+1, qcd_mc_est.at(i));
            hp_mm_qcd->SetBinError(i+1, qcd_mc_est_err.at(i));
        }        
    }

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    // creating the plot
    c = new TCanvas(plot_title, plot_title,1);
    c->SetTickx();
    c->SetTicky();

    if (logy) c->SetLogy();

    hp_mc_qcd->Divide(hp_mm_qcd); set_jet_binlabel(hp_mc_qcd, true);

    hp_mc_qcd->GetXaxis()->SetTitle("number of jets");
    hp_mc_qcd->GetYaxis()->SetTitle("MC QCD / MC MM");
    hp_mc_qcd->GetXaxis()->SetLabelSize(.03);
    hp_mc_qcd->GetYaxis()->SetLabelSize(.03);
    hp_mc_qcd->GetXaxis()->SetTitleSize(.03);
    hp_mc_qcd->GetYaxis()->SetTitleSize(.03);
    hp_mc_qcd->GetXaxis()->SetTitleOffset(1.5);
    hp_mc_qcd->GetYaxis()->SetTitleOffset(1.5);
    hp_mc_qcd->SetMarkerStyle(20);

    hp_mc_qcd->Draw();

    /*
    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_mc_qcd_ratio,"MC_QCD / MC_MM","p");
    legend->Draw();
    */

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}

TCanvas * cal_ctrl_syserr(vector<TString> hist_fpaths, Double_t d0sig_norm, const int ijet) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    Double_t qcd_norm = 0;

    const Int_t n = hist_fpaths.size();

    Double_t qcd[n];
    Double_t qcd_err[n];
    Double_t d0sig_err[n];
    Double_t d0sig[n];

    ostringstream stm_title;
    stm_title << "sys_err_rctl_j" << ijet;
    if ( hist_fpaths[0].Contains("_btag") ) stm_title << "_btag";
    if ( hist_fpaths[0].Contains("_tria") ) stm_title << "_tria";
    TString plot_title = TString(stm_title.str());

    for (int i=0; i<n; i++) {

        TFile * f     = new TFile( hist_fpaths[i] );
        TH1D  * h_qcd = (TH1D*) f->Get("data10_njets_qcd_estat_j0"); h_qcd->SetDirectory(0);
        TAxis * x     = h_qcd->GetXaxis();

        Double_t qcd_4j     = 0.;
        Double_t qcd_4j_err = 0.;
        for(int j=1; j <= h_qcd->GetNbinsX(); j++) {
            Double_t N_jets    = x->GetBinLowEdge(j);
            Double_t N_qcd     = h_qcd->GetBinContent(j);
            Double_t N_qcd_err = h_qcd->GetBinError(j);

            //cout << N_jets << " " << N_qcd << " " << N_qcd_err << endl;

            switch(ijet) {
                case 0:
                    qcd_4j     += N_qcd;
                    qcd_4j_err += TMath::Power(N_qcd_err,2);
                    break;

                case 4:
                    if (N_jets >= 4) {
                        qcd_4j     += N_qcd;
                        qcd_4j_err += TMath::Power(N_qcd_err,2);
                    }
                    break;

                default:
                    if (N_jets == ijet) {
                        qcd_4j     += N_qcd;
                        qcd_4j_err += TMath::Power(N_qcd_err,2);
                    }
                    break;
            }
        }

        qcd_4j_err = TMath::Sqrt(qcd_4j_err);

        qcd[i]       = qcd_4j;
        qcd_err[i]   = qcd_4j_err;
        d0sig[i]     = parse_d0sig( hist_fpaths[i] );
        d0sig_err[i] = 0.5;

        if ( d0sig[i] == d0sig_norm ) {
            qcd_norm = qcd[i];
        }

        f->Close();

        delete f;
    }

    // printing the variations of QCD estimation vs. control region
    Double_t stdev = 0.;
    cout << "========================" << endl;
    for( int i=0; i<n; i++) {

        stdev += TMath::Power( qcd[i] - qcd_norm, 2);

        Double_t sys_err = TMath::Abs( qcd[i] - qcd_norm ) / qcd_norm;
        cout << d0sig[i] << " " << qcd[i] << " " << qcd_err[i] << " " << sys_err << endl;
    }
    stdev = TMath::Sqrt(stdev);
    Double_t stdev_ratio = stdev * 100.0 / qcd_norm;
    cout << "========================"                           << endl;
    cout << "std. dev: " << stdev << " (" << stdev_ratio << "%)" << endl;
    cout << "========================"                           << endl;

    // creating the plot
    c = new TCanvas("control",plot_title,200,10,700,500);
    //c->SetFillColor(42);
    //c->SetGrid();
    c->GetFrame()->SetFillColor(21);
    c->GetFrame()->SetBorderSize(12);

    // create the TGraphErrors and draw it
    g1 = new TGraphErrors(n,d0sig,qcd,d0sig_err,qcd_err);
    g1->SetTitle("statistical uncertainty");
    g1->GetXaxis()->SetTitle("d0 significance");
    g1->GetYaxis()->SetTitle("number of events");
    g1->GetXaxis()->SetLabelSize(.03);
    g1->GetYaxis()->SetLabelSize(.03);
    g1->GetXaxis()->SetTitleSize(.03);
    g1->GetYaxis()->SetTitleSize(.03);
    g1->GetXaxis()->SetTitleOffset(1.5);
    g1->GetYaxis()->SetTitleOffset(1.5);
    g1->SetFillColor(42);
    g1->Draw("AE3");

    g2 = new TGraphErrors(n,d0sig,qcd,d0sig_err,qcd_err);
    g2->SetTitle("qcd estimation 4 jets incl.");
    g2->SetMarkerColor(kBlack);
    g2->SetMarkerStyle(20);
    g2->Draw("SAME PZ");

    // drawing horizontal lines showing the deviations based on the norminal value
    TLine *l1 =  new TLine(d0sig[0], qcd_norm+stdev, d0sig[n-1], qcd_norm+stdev);
    TLine *l2 =  new TLine(d0sig[0], qcd_norm-stdev, d0sig[n-1], qcd_norm-stdev);

    l1->SetLineWidth(2); l2->SetLineWidth(2);
    l1->SetLineColor(2); l2->SetLineColor(2);
    l1->SetLineStyle(2); l2->SetLineStyle(2);
    l1->Draw("SAME");    l2->Draw("SAME");

    ostringstream stm;
    stm << "#sigma_{d0sig=" << d0sig_norm << "} = " << stdev;
    TString stdev_str = TString( stm.str() );
    myText(0.65,0.8, kRed, stdev_str);

    c->Update();

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    return c;
}

void cal_mc2data_sfactor(TString hist_fpath) {

    TFile *f = new TFile(hist_fpath);

    TH1D *hqcd_data = f->Get("data10_njets_qcd_estat_j0"); hqcd_data->SetDirectory(0);
    TH1D *hqcd_mc   = f->Get("mc_qcd_njets_j0")          ; hqcd_mc->SetDirectory(0);
    TH1D *hqcdbb    = f->Get("mc_qcdbb_njets_j0")        ; hqcdbb->SetDirectory(0);
   
    hqcd_mc->Add(hqcdbb);
   
    vector<Double_t> qcd_data_est;  qcd_data_est.resize(5);
    vector<Double_t> qcd_mc_est;    qcd_mc_est.resize(5);
    vector<Double_t> qcd_est_err;   qcd_est_err.resize(5);
   
    qcd_est_by_jets(hqcd_mc,   &qcd_mc_est,   &qcd_est_err);
    qcd_est_by_jets(hqcd_data, &qcd_data_est, &qcd_est_err);

    const int s = (int) qcd_data_est.size();

    for (int i=0; i<s; i++) {
        Double_t r = qcd_data_est.at(i) / qcd_mc_est.at(i);
        cout << "ijet=" << i << " " << qcd_data_est.at(i) << " " << qcd_mc_est.at(i) << " " << r << endl;
    }
}

void cal_contamination_syserr(TString hist_fpath) {

    TFile *f = new TFile( hist_fpath );

    TH1D *h_qcd_n       = (TH1D*) f->Get("data10_njets_qcd_estat_j0"); h_qcd_n->SetDirectory(0);

    TH1D *h_njets_tight = (TH1D*) f->Get("data10_njets_tight_j0"); h_njets_tight->SetDirectory(0);
    TH1D *h_njets_loose = (TH1D*) f->Get("data10_njets_loose_j0"); h_njets_loose->SetDirectory(0);

    TH1D *h_my_njets_tight = h_njets_tight->Clone();
    TH1D *h_my_njets_loose = h_njets_loose->Clone();

    h_my_njets_tight->Reset();
    h_my_njets_tight->Add(h_njets_tight);

    h_my_njets_loose->Reset();
    h_my_njets_loose->Add(h_njets_loose);

    // recalcuate fake efficiency and make new QCD estimation
    vector<Double_t> rec_eff_fake_p = eff_fake_contaminate( hist_fpath, 0, '+' );
    TH1D *h_qcd_p = PlotQCD(h_my_njets_loose, h_my_njets_tight, rec_eff_fake_p[0], rec_eff_fake_p[1], 4, 0, 0);

    h_my_njets_tight->Reset();
    h_my_njets_tight->Add(h_njets_tight);

    h_my_njets_loose->Reset();
    h_my_njets_loose->Add(h_njets_loose);

    // recalcuate fake efficiency and make new QCD estimation
    vector<Double_t> rec_eff_fake_m = eff_fake_contaminate( hist_fpath, 0, '-' );
    TH1D *h_qcd_m = PlotQCD(h_my_njets_loose, h_my_njets_tight, rec_eff_fake_m[0], rec_eff_fake_m[1], 4, 0, 0);

    // count QCD estimation by jets
    vector<Double_t> qcd_n_est;    qcd_n_est.resize(5);
    vector<Double_t> qcd_p_est;    qcd_p_est.resize(5);
    vector<Double_t> qcd_m_est;    qcd_m_est.resize(5);
    vector<Double_t> qcd_est_err;  qcd_est_err.resize(5);
   
    qcd_est_by_jets(h_qcd_n, &qcd_n_est, &qcd_est_err); 
    qcd_est_by_jets(h_qcd_p, &qcd_p_est, &qcd_est_err);
    qcd_est_by_jets(h_qcd_m, &qcd_m_est, &qcd_est_err);

    const int s = (int) qcd_p_est.size();

    Double_t r_p;
    Double_t r_m;
    for (int i=0; i<s; i++) {
        r_p = TMath::Abs(qcd_p_est.at(i) - qcd_n_est.at(i)) / qcd_n_est.at(i);
        r_m = TMath::Abs(qcd_m_est.at(i) - qcd_n_est.at(i)) / qcd_n_est.at(i);
        cout << "ijet=" << i << " " << qcd_n_est.at(i) << " " << qcd_m_est.at(i) << "(" << r_m << ") " << qcd_p_est.at(i) << "(" << r_p << ")" << endl;
    }
}
